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ON THE ECCENTRICITY DISTRIBUTION OF EXOPLANETS FROM RADIAL VELOCITY SURVEYS 
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ABSTRACT 

We investigate the estimation of orbital parameters by least-x 2 Keplerian fits to radial velocity (RV) data us- 
ing synthetic data sets. We find that while the fitted period is fairly accurate, the best-fit eccentricity and M p sin i 
are systematically biased upward from the true values for low signal-to-noise ratio K/cr < 3 and moderate 
number of observations N b s < 60, leading to a suppression of the number of nearly circular orbits. Assuming 
intrinsic distributions of orbital parameters, we generate a large number of mock RV data sets and study the 
selection effect on the eccentricity distribution. We find the overall detection efficiency only mildly decreases 
with eccentricity. This is because although high eccentricity orbits are more difficult to sample, they also have 
larger RV amplitudes for fixed planet mass and orbital semi-major axis. Thus the primary source of uncer- 
tainties in the eccentricity distribution comes from biases in Keplerian fits to detections with low-amplitude 
and/or small A/ b s , rather than from selection effects. Our results suggest that the abundance of low-eccentricity 
exoplanets may be underestimated in the current sample and we urge caution in interpreting the eccentricity 
distributions of low-amplitude detections in future RV samples. 
Subject headings: planetary systems - techniques: radial velocities 



1. INTRODUCTION 

With the rapid increase in the rate of exoplanet detections, it 
has become feasible to study their statistical properties (for a 
recent review, see Udry & Santos 2007). The distributions of 
their orbital parameters and correlations with host star prop- 
erties are crucial for our understanding of planet formation. 
Up to Feb 2008, over 200 exoplanets have been announced, 
most of which were detected by the radial velocity (RV) tech- 
nique. Among the parameters that can be derived from RV 
data, the orbital eccentricity has a somewhat unexpected dis- 
tribution, with an extended tail of high eccentricities (<? > 0. 1). 
Although the exact form of this eccentricity distribution is still 
uncertain to some extent, especially at the lowest eccentrici- 
ties (i.e., comparing the distribution in Butler et al. 2006 and 
in the most recent catalog), it is clear that this distribution is 
quite different from that of the Solar system. There have al- 
ready been several theoretical attempts to explain such an ec- 
centricity distribution (e.g., Tremaine & Zakamska 2004 and 
references therein; Juric & Tremaine 2007; Ford et al. 2007; 
Zhou et al. 2007). 

However, the largest current exoplanet catalog (e.g., Bul- 
ter et al. 2006 and their updates, hereafter Bulter06) is by no 
means homogeneous. Survey strategies, selection biases in 
the RV technique and uncertainties in best-fit orbital solutions 
could all bias the intrinsic distributions of orbital parameters, 
especially when taken together. Eccentric orbits have larger 
amplitudes than circular orbits when holding other parameters 
fixed, while failure to (time) resolve the perihelion approach 
can lead to non-detection for an eccentric orbit system. These 
two effects are believed to roughly cancel but detailed simu- 
lation is needed (e.g., Endl et al. 2002; Cumming 2004). It 
is also well known that the errors in least-x 2 Keplerian so- 
lutions become asymmetric for noisy data (e.g., Ford 2005; 
Butler et al. 2006), which is especially true for eccentricity. 
On the other hand, for circular orbits the eccentricity can only 
be scattered upward by error, so an obvious bias exists for 
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circular orbits. As the RV surveys progress, more and more 
low amplitude (low signal-to-noise ratio) detections will be 
reported. Significant uncertainties remain in the distributions 
of orbital parameters for those exoplanets using best-fit or- 
bital solutions. There are already some cases where the orbit 
eccentricity is poorly constrained by the RV data (e.g., Jones 
et al. 2006). 

The purpose of the paper is to explore processes that might 
distort measurements of the shape of the intrinsic eccentricity 
distribution using synthetic data. First we investigate the relia- 
bility of single Keplerian fits to mock RV data sets as function 
of signal-to-noise ratio, and see if any bias arises from Kep- 
lerian fitting to noisy data (©. Second, we construct statis- 
tical models of orbital parameter distributions and pass them 
through a simulated detection pipeline, in order to see if there 
is any serious selection effect on eccentricity ($3). We discuss 
our results and apply them to the current RV planet sample in 
<2] and we summarize our findings in $5] 



2. KEPLERIAN FITTING TO MOCK RV DATA 

We consider single planet detections throughout the pa- 
per. Assuming a true orbit (period P, eccentricity <?, minimal 
planet mass M p sin i, argument of periastron u, and tempo- 
ral offset fo) and the host star mass M», we generate mock RV 
data sets with Gaussian errors a as perturbations. For simplic- 
ity <7 is assumed to be the same for all data points in a set, but 
we relax this assumption later (§|4). In practice, a comes from 
RV measurement errors, stellar jitter and possible existence 
of multiple planets (e.g., Ford 2006) 1 . Following Cumming 
(2004), we define the signal-to-noise ratio 1Z = K/a where K 
is the semi-amplitude of the theoretical RV curve, and in the 

1 We note that in general the additional reflex motion induced by unknown 
secondary planets cannot be treated as independent perturbations. Second, 
for simplicity we do not consider secular trends in the orbital motion of the 
primary planet induced by unseen companions, which are typically on much 
longer timescales than the period and observation time baseline considered 
here. 
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FIG. 1. — Effects of signal-to-noise ratio and number of observations on the best-fit orbital parameters. Filled circles are median values; the true orbital 
parameters are denoted as horizontal lines (see the text for details), where u and temporal offset tQ are chosen at random. Error bars show standard deviations. 
From left to right, N obs = 10, 20, 60. 

to this point in ^3] The mock observational time span T was 
set to be two true periods, but we find almost identical results 
as long as the period is shorter than the time span. Finally, we 
simulated three values of Af bs = 10, 20, 60. 

Not all attempts to find a Keplerian fit will result in a de- 
tection. In an actual radial velocity search for planets, some 
threshold on the false alarm probability (FAP), the probability 
that a signal power will arise purely from noise, of a Keple- 
rian fit, is normally used. Throughout the paper we adopt a 
fiducial FAP = 10~ 2 , which is essentially always satisfied for 
large K/a and A^,b s . There are several ways to estimate the 
FAP using Monte Carlo simulations or analytical formulae. 
For simplicity we have followed the analytical procedure de- 
scribed in Cumming (2004) to estimate the FAP: Given the 



Here the amplitude K is defined in terms of the intrinsic prop- 
erties of a planet's orbit, since we "know" the true orbital pa- 
rameters in advance. 

In addition to the signal-to-noise ratio 1Z, another important 
quantity is the number of observations A^,b s - Clearly the larger 
1Z and A^,bs are, the better the quality of the data set. 

We consider planets where the orbital period is shorter than 
the observational time span. Orbits with periods longer than 
the time baseline are difficult to detect due to limited sampling 
and lower amplitude (e.g., Cumming 2004). Even if detected, 
the best-fit orbital solution would often be very different from 
the true one (Shen et al. 2008, in preparation). The obser- 
vation times are randomly generated within the observational 
time span, which somewhat mimics the effects of realistic ra- 
dial velocity sampling. 

In Fig. Q]we plot the relations of fitted eccentricity, period, 
minimal planet mass and semi-amplitude to their true values, 
as functions of 1Z, for a specific case with M* = 1.1 M Q , P = 50 
days, M p sin i = 1 My, and eccentricity values e = 0.0-0.5 at 
intervals of Ae = 0.1. Each set of true orbits was used to 
generate 5000 mock RV data sets with randomly chosen u> and 
temporal offset fo, and with Gaussian noise determined from 
the assumed signal-to-noise ratio. We used the Levenberg- 
Marquardt method (e.g., Press et al. 1992) to minimize x 2 > 
and the true orbital parameters were used as an initial guess 
for the x 2 fitting. Using the true orbital parameters as the 
initial guess greatly speeds up the convergence and increases 
the rate of successful fits. In the realistic case when we do not 
know the actual orbital parameters, the success of a fit also 
depends on an appropriate initial guess. We will come back 



least-x value xk ep fr° m the Keplerian fit, and the least-x 
value Xmean from the fit of a constant to the data, we cal- 
culate a power zo using eqn. (7) in Cumming (2004). An 
estimate of the FAP is then eqn. (5) in Cumming (2004): 
FAP = 1 - [ 1 - Prob(z >zo)] M , where M^T/ T min is the num- 
ber of independent frequencies with T the observational time 
span and T m \ n = 2 days the lower bound of period searched in 
the Keplerian fitting, and Prob(z > zo) is the probability dis- 
tribution given by eqn. (8) in that paper. 

All three panels in Fig. Q]show that the fitted period is fairly 
accurate. However, for eccentricity and M p sin / (and also the 
amplitude K following eqn. [TJ the median fitted values are 
biased upward, and this behavior remains for true eccentric- 
ities as large as eo ~ 0.8 in our simulations. This is not a 
failure of the \ 2 minimization method. In Fig. [2] we show 
two examples of Keplerian fits where both RV data sets prefer 
a more eccentric orbit. The biases in e and M p sini, as well 
as the scatter around the median, decrease when 7Z or iV bs 
increases, as would be naively expected. 

These fitting biases result from the fact that, in general, the 
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FIG. 2. — Left: An example of Keplerian fit to mock data with K/cr = 5 and N B b s ■■ 
original circular orbit; but the fitted period and M p sinr are fairly accurate. Right. An example of large eccentricity orbit fits with K/cr = 2 and N bs = 20 where 
eo = 0.5. In both cases the dashed lines are the RV curves expected from the true orbital parameters while the solid lines are best-fit Keplerian orbits to the data. 

least-x 2 solution is not an unbiased estimator of the true pa- 



rameters for the non-linear Keplerian model even if the data 
errors are Gaussian. In particular, the probability distribution 
V(6\6) is not Gaussian, where 9 is the set of true parameter 
and 9 is the best-fit orbital solution; and it only approaches a 
Gaussian distribution when the errors are small. 

To fully take account these biases we need to know 7 5 (f?|f?) 
for different combinations of 9, K/a and iVobs, which can be 
explored numerically as in Fig. [Tj However, for low intrin- 
sic eccentricities (eo < 0.5), the apparent fitting bias in ec- 
centricity can be quantitatively understood as follows. While 
the error in the fitted eccentricity a e is not Gaussian when 
a e is comparable to e, the fitted x = ecosw and y = esinu 
preserve Gaussianity much better provided that the intrinsic 
eccentricity is small (e.g., Butler et al. 2006). Assuming 
their fitted values are Gaussian distributions around true val- 
ues xrj = erjcoswo and yo = eosinojo with dispersion ao, the 
probability distribution of the fitted eccentricity e is given by 



V(e)(x / exp 



(x-x ) (y-yo) 



2a% 



5(e- \J x 2 +y 2 )dxdy 



oceexp 



2a\ 



I eeo_ 

u 



(2) 



where Iq is the modified Bessel function of the first kind. For 
circular orbits (e ( > = 0) it is a Rayleigh distribution; for ec- 
centric orbits and at the low-dispersion limit (ctq 0) it is 
approximately a Gaussian distribution. 

Fig. [3] shows an example of [P,eo,M p sin;,a;o,fo,M»] = 
[50,0.1, l,7r/4, 10, 1.1] with JV obs = 60. For a high signal-to- 
noise ratio, K/a = 5 (lower panels), both distributions of fitted 
eccentricity and ecosw (esinw) are approximately Gaussian. 
For a low signal-to-noise ratio, K/a = 2 (upper panels), the 
eccentricity distribution significantly deviates from a Gaus- 
sian while ecosw and e sin u are still distributed in an approx- 
imately Gaussian form. We plot the predictions of fitted ec- 
centricity distribution from eqn. © as solid lines, using the 
best-fit Gaussian dispersion of the e cos ui/e sin ui distributions. 
In this way we can also understand the upward median biases 
even at intrinsic eccentricities as high as ~ 0.8, although the 
fitted ecoscj and esinuj are no longer Gaussian at such high 
eccentricities-our simulations themselves are more revealing 
in this case. 

A conservative empirical criterion (based on our numer- 
ical experiments) for a median bias in fitted eccentricity 
((e) - eo) ^ 0.05 is given by the following joint constraint on 



.Q 

O 



0.5 
0.4 

0.3 

0.2 

0.1 

0.0 
0.5 

0.4 
0.3 



K/<7 = 2 



N„ 



= 60 





o 0.2 

0.1 

0.0 



K/<7 = 5 N 



= 60 



J 




<e> = 0.1 10 



-0.5 0.0 0.5 
ecoscj.esincj 



0.0 0.1 0.2 0.3 0.4 0.5 0.6 
Eccentricity 

FIG. 3. — Examples of fitted distributions becoming non-Gaussian when 
the signal-to-noise ratio decreases (see text for the intrinsic parameters). Left: 
Distributions of fitted e cos u> (solid histogram) and e sin uj (dotted histogram). 
Right: Distributions of fitted eccentricities. The median value of fitted eccen- 
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from eqn. {2) and the simulated distributions is due to the fact that the distri- 
butions of ecoscj and esinu are not precisely Gaussian. 



K/a and A^oi 



> 15 



(3) 



which works reasonably well for 10 < A^,bs < 200. This is also 
an approximate criterion that the median bias in M p sin i is less 
than 10%. Note again here that K is defined using intrinsic 
planet properties. In practice, K can only be estimated from 
the best-fit solution and will on average be overestimated a 
bit at the low-amplitude end, according to the fitting biases 
discussed here (e.g., Fig. Q]i. If using Kf& instead, we found a 
good approximation is to replace equation (O with 



(KmMnH* > 17 , 



(4) 



which compensates the fitting bias of amplitude at the low 
K/a end. 



3. 



STATISTICAL MODELS OF ORBITAL PARAMETER 
DISTRIBUTIONS 
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FIG. 4. — Colored lines show the detection efficiencies as a function of 
eccentricity within a given bin of Af b s ■ The long-dashed cyan line show the 
detection efficiency for fixed signal-to-noise ratio. The solid and dashed his- 
tograms show the true eccentricity distributions for all the mock orbits and 
the detected orbits, with the same peak normalization. 

We now investigate the selection bias on eccentricity us- 
ing synthetic RV data sets generated from model distributions 
of orbital parameters. The detailed comparison of different 
model distributions and constraints from observational data 
will be presented in a paper in preparation (Shen et al. 2008). 
Our purpose here is to investigate the possible selection bias 
on eccentricity, therefore essentially any well behaved model 
distributions will do. Nevertheless, we choose a model with 
physically plausible distributions, as we describe below. 

The period distribution is a log-normal with log mean of 10 
yrs and dispersion of 1 dex. The planet mass M p distribution 
is a power-law ndM p oc M" 1 within [0. 1 , 100] Mj; this extends 
well into the "brown dwarf desert" due to the upper limit of 
100 Mj. However such objects are too rare to have any sta- 
tistical impact on our results. The orientation of the planet 
orbital plane is random, as is the argument of periastron u> 
within [0,27r] and the temporal offset to within one period. 
For the eccentricity distribution we chose a model that peaks 
at e = and diminishes to zero at e = 1 : 



V(e)de oc 



1 



(l + e) a 2 a 



de 



(5) 



and we assume a = 4 for this specific study. Our choice of 
this particular model distribution is because the peak of the 
observed eccentricity distribution has shifted from e w 0.2 to 
e w using the latest exoplanet catalog 2 . This model distri- 
bution of eccentricity is shown as the solid histogram in Fig. 
|4] Finally we specify an observational time baseline of T = 10 
years and an overall radial velocity error a = 10 ms" 1 for each 
simulated observation; the number of observations of each 
object is randomly distributed from A^bs = 20 to A^bs = 100. 
These parameters are intended to be typical of currently RV 
surveys (e.g., Cumming et al. 2008). 

We generate 100,000 mock orbits of the model distributions 
and pass them through the Keplerian detection pipeline with 
a detection threshold FAP = 0.01. For the initial guess of 

2 http://www.exoplanets.org 
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FIG. 5. — Low K/a detections for the statistical model (8157 planets de- 
tected). Higher and lower histograms show the distributions of the intrinsic 
and best-fit eccentricities of these detected planets. Shown here is the bias 
caused by Keplerian fitting: ~ 40% of the nearly-circular (e < 0.1) orbits are 
erroneously assigned to more eccentric orbits by their best fits. 



orbital solution, we again take advantage of our mock data 
sets by using the true orbital parameters. In the realistic case, 
the initial guess for period can be identified using traditional 
Lomb-Scargle periodogram (Lomb 1976; Scargle 1982), but 
appropriate initial guesses for [e, oj, to] are crucial for the con- 
vergence of Keplerian fits for orbits with high eccentricities, 
where the fit can easily be trapped in a local minimum due 
to the complicated \ 2 space or fail to converge at all, and 
hence the planet is undetected. We have run a substantial 
number of tests, where we tried a fine grid of [e,w,fo] as ini- 
tial guesses in each Keplerian fit. While time consuming, we 
found that in the vast majority cases the Keplerian fits find 
the same minimum as using the true parameters as the initial 
guess. Serendipitous aliasing cases do exist, where the global 
minimum solution deviates catastrophically from the true or- 
bital parameters, and which are for the most noisy data sets. 
Hence in the realistic case, special care must be taken to find 
the global minimum efficiently. 

In Fig. 2]we show the overall detection efficiency (the frac- 
tion of detected orbits) as function of eccentricity in filled cir- 
cles, and different line types represent the detection efficiency 
when binned in N b s - Clearly increasing Aob s increases the de- 
tection efficiency. However, the overall detection efficiency 
only slowly decreases with eccentricity. This is because al- 
though high eccentricity orbits are more difficult to sample, 
they also have larger amplitude K on average. When we plot 
the detection efficiency with constant amplitude Kj a ps 2, we 
see a decrease for e > 0.6, as shown by the cyan long-dashed 
line in Fig. HI which is consistent with previous studies (e.g., 
Endl et al. 2002; Cumming 2004). Therefore, there is no 
strong selection bias against eccentricity, as demonstrated by 
the similar (true) eccentricity distributions of the mock orbits 
and the detected ones in Fig. |4] Note, however, the absolute 
value of the overall detection efficiency here depends on the 
model distributions and the quality of the survey (noise level 
and number of observations). 

However, as we have shown in ^2j \ 2 Keplerian fits will 
have non-negligible effects for low amplitude detections and 
therefore bias the underlying eccentricity distribution. Al- 
though the detection efficiency decreases rapidly when K/a 
decreases, there are also many more objects in the low am- 
plitude regime given our model distributions and noise level. 
We show this bias effect in Fig. |5j where the blue and red 
histograms show the distributions of true and fitted eccentric- 
ities for low amplitude detections with K/a < 3. The small 
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FIG. 6. — The effects of Keplerian fitting on the eccentricity distributions. The solid and dashed histograms are the distributions of intrinsic and fitted 
eccentricities in different K/a and N \, s bins. All planets are detected by the Keplerian fitting pipeline, and the numbers of detected planets are listed at the 
top-right corner of each panel. Those bins with the lowest K/cr and N ^ s are the most seriously affected by the fitting bias. 



pile up of fitted e > 0.9 orbits is caused by detections with 
period longer than the observation time where the Keplerian 
fits failed badly (Shen et al. 2008, in preparation). Other than 
that, the overall effect of the Keplerian fitting bias is the trans- 
fer of true low-eccentricity orbits to fitted high-eccentricity 
ones, therefore distorting the underlying eccentricity distri- 
bution. As expected from our simulations in §|2] this effect 
diminishes when we shift to high amplitude and/or large A^bs 
detections, as shown in Fig. [6] 

4. DISCUSSION 

Our results lead naturally to the question of how much bias 
in the eccentricity distribution is to be expected in the Butler et 
al. (2006) catalog. To test this we take their updated catalog 
as of January 26, 2008. We exclude those planets with or- 
bital periods less than 20 days (i.e., those probably have been 
tidally circularized). Over half of the planets in that catalog 
do not have published RV data and measurement errors/stellar 
jitter information, so we use the RMS scatter of the model fit 
as an estimate of the noise a. We also take the amplitude 
K as the best-fit value Ret. There are 162 planets with K^, 
RMS scatter and A^,b s available. In Fig. [7] we show the distri- 
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FIG. 7. — Inspection of the exoplanet catalog as of January 26, 2008 
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the fitting bias (e.g., eqn. |4}. 

bution of (Kst/a)N^ defined in equation for these plan- 
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ets. There are 19 planets with (K&t/o^NgL < 17. Thus sta- 
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8 . — An example of the Kepler fits to RV data in the Butler06 catalog 
14729b). This particular case has iV bs = 42 and median Kfn/cr ks 3.5, 
a is estimated from the published RV measurement errors and stellar 



Fig. 
(HD 1 
where 
jitter. 

tistically speaking, only ~ 10% of the planets are affected in 
the current sample, which are predominately low-amplitude 
detections. 

To test the effects of unequal errors and realistic sampling 
we take HD 1 14729b in the Butler06 catalog as an example. It 
has a best-fit eccentricity 0.167 ±0.055 with Kfn/a « 3.5 and 
/V bs = 42 where we have estimated a using their published RV 
measurement errors and stellar jitter. The actual RV data and 
their best-fit solution is shown in Fig. [8] Using the true time 
series and data uncertainties derived from measurement errors 
and stellar jitter for this system, but with simulated radial ve- 
locities, we find that the probability of a fitted eccentric orbit 
with e > 0.15 arising from a circular orbit is ~ 8%, and in- 
creases to 20% for a fitted e > 0.1, consistent with our results 
in ^2]with constant uncertainties and random sampling. 

5. CONCLUSIONS 

We have performed simulations of planet detection and or- 
bital parameter fits using mock radial velocity data sets. Two 
effects that may affect the intrinsic eccentricity distribution 
are considered: selection bias on eccentricity, and fitting bias 
in the best-fit orbital solution. 

We find that selection bias on eccentricity is negligible, as 
long as the fitting routine is efficient in finding the global so- 
lution. In a realistic survey, this requires a thorough search 
in the parameter space, and/or with advanced algorithms op- 
timized for the search. Our finding is not in conflict with pre- 
vious studies (e.g., Endl et al. 2002; Cumming 2004), which 



claim that the detection efficiency decreases for e > 0.6. This 
is because in previous studies, the detection efficiency as a 
function of eccentricity is estimated at fixed K while in our 
study K is larger for more eccentric orbits when other param- 
eters are fixed. On the other hand, we find that for detections 
with low signal-to-noise ratio and small number of observa- 
tions, the best-fit eccentricity is biased upward in the median 
value, which then gives rise to a change in the eccentricity 
distribution from the intrinsic one. 

Inspection of the current exoplanet sample shows only ~ 
10% are likely to be significantly affected by the Keplerian 
fitting bias. However, future radial velocity surveys will con- 
tain an increasing number of low amplitude detections if the 
number of low mass, large semi-major axis exoplanets grows 
as rapidly as suggested by extrapolation of current results sug- 
gest (e.g., Udry & Santos 2007). When the sample of low am- 
plitude detections (on average less massive planets and more 
circular orbits) is large enough for statistical study, the bias in 
the best-fit eccentricity described here must be taken into ac- 
count. Also, in individual cases where an accurate estimate of 
eccentricity is required, i.e., in modeling the habitable zone or 
tidal heating issues, either high quality RV data or constraints 
from other observations are required to support reliable con- 
clusions. In the mean time, there is also a need for a more 
statistically sophisticated understanding of the uncertainties 
of derived orbital solutions (e.g., Ford 2005; 2006). 

Our results suggest that the intrinsic eccentricity distribu- 
tion may be even more peaked at e m than the current 
observed distribution. Some planet-planet scattering models 
tend to produce a Rayleigh-distribution of eccentricities (e.g., 
Juric & Tremaine 2007; Ford & Rasio 2007) with reduced 
circular orbits, therefore certain eccentricity damping mecha- 
nism such as interactions with a protoplanetary disk may be 
required to reconcile these models with observations. 
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